# Project Info----------------------------------------------------------------------
#  Project: Outgroup Avoidance
#  Purpose: Analyze Facebook Field Experiment 
#  Outputs: Figure 6, Tables S2-5, Figure S4

# load all relevant packages--------------------
library("tidyverse")   
library("stargazer")   
library("texreg")      
library("estimatr")    
library("stdidx")   # Standardizing variables
library("modelsummary") # Model summaries & tables
library("kableExtra")   # Formatting tables


# load data-------
stdy1_dta <- read.csv("study1.csv")


# Main Text------
## Figure 6-----
therm_result <-
  lm_robust(stdize(therm_arabs_post_num) ~
              pooled + arab_therm_pre + arab_therm_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Feeling \nThermometer",
         group = "Attitudes") %>% 
  filter(.,
         term == "pooled")

empathy_result <- 
  lm_robust(stdize(empathy_index_post) ~
              pooled + empathy_index_pre + empathy_index_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Empathy \nIndex",
         group = "Attitudes") %>% 
  filter(.,
         term == "pooled")

peace_result <-
  lm_robust(stdize(arabs_peace) ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Peace",
         group = "Attitudes") %>% 
  filter(.,
         term == "pooled")

exposure_attitudes_result <-
  lm_robust(stdize(exposure_index_post) ~
              pooled + exposure_positive_pre +exposure_positive_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Exposure \nIndex",
         group = "Attitudes") %>% 
  filter(.,
         term == "pooled")

video_result <-
  lm_robust(stdize(video) ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Partner \nOrganization \nVideo",
         group = "Behaviors") %>% 
  filter(.,
         term == "pooled")

follow_result <-
  lm_robust(stdize(page_like) ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Follow \nPartner \nOrganization",
         group = "Behaviors") %>% 
  filter(.,
         term == "pooled")

exposure_bhvr_result <-
  lm_robust(stdize(bhvr_indx) ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Exposure \nBehavior \nCombined",
         group = "Behaviors") %>% 
  filter(.,
         term == "pooled")

##Plot results------
main_models <-
  rbind(therm_result, empathy_result, peace_result,
        exposure_attitudes_result, exposure_bhvr_result,
        follow_result, video_result) %>% 
  mutate(.,
         b = paste0("b=", formatC(estimate, format = "f", digits = 3)),
         p = paste0("p=", formatC(p.value, format = "f", digits = 3)),
         bp = paste0(b, ", ", p))

ggplot(main_models, aes(x=estimate, y=outcome)) + 
  geom_point(color = "dodgerblue4")+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1, color = "dodgerblue4") +
  facet_wrap(~group,
             scales = "free") +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             linewidth = 0.4) +
  labs(x = "",
       y = "")+
  geom_text(aes(label= bp), 
            colour = "dodgerblue4",
            position = position_dodge(width = .6),
            hjust = .4,
            vjust = -0.6,
            size = 4,
            show.legend = FALSE) +
  scale_color_manual(values = c("dodgerblue4"))+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 16),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())




#Appendix-----
##Table S2------
atrt_pre <-
  lm_robust(no_pre ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)

atrt_post <-
  lm_robust(no_post ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)

atrt_pre_post <-
  lm_robust(no_prepost ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)

atrt_therm <-
  lm_robust(no_post_therm ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)

atrt_emp <-
  lm_robust(no_post_empathy ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)

atrt_peace <-
  lm_robust(no_post_peace ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)

atrt_post_expose <-
  lm_robust(no_post_expose ~
              pooled,
            fixed_effects = block_id,
            data = stdy1_dta)


models <- 
  list("Pre" = atrt_pre, 
       "Post" = atrt_post, 
       "Pre-Post" = atrt_pre_post, 
       "Therm" = atrt_therm, 
       "Empathy" = atrt_emp,
       "Peace" = atrt_peace,
       "Exposure" = atrt_post_expose)

modelsummary(models,
               output ="latex",
               title = "Correlation of Treatment with Non-Response to Surveys and Specific Items\\label{attrition}",
               coef_map = c(pooled = "Treatment"),
               gof_map = c("nobs", "r.squared"),
               stars = c("*"=0.05))  %>%
  add_header_above(c("Outcome: Non Response = 1, Response = 1" = 8))



##Table S3------
cov_post <-
  lm_robust(no_post~
              male + age_num + religiosity_num + therm_arabs_num,
            data = stdy1_dta)

cov_therm <-
  lm_robust(no_post_therm~
              male + age_num + religiosity_num + therm_arabs_num,
            data = stdy1_dta)


cov_emp <-
  lm_robust(no_post_empathy~
              male + age_num + religiosity_num + therm_arabs_num,
            data = stdy1_dta)


cov_peace <-
  lm_robust(no_post_peace~
              male + age_num + religiosity_num + therm_arabs_num,
            data = stdy1_dta)

cov_exp <-
  lm_robust(no_post_expose ~
              male + age_num + religiosity_num + therm_arabs_num,
            fixed_effects = block_id,
            data = stdy1_dta)



cov_models <- 
  list("Post" = cov_post, 
       "Therm" = cov_therm, 
       "Empathy" = cov_emp,
       "Peace" = cov_peace,
       "Exposure" = cov_exp)

modelsummary(cov_models,
               output ="latex",
               title = "Correlation of Covariates with Non-Response to Surveys and Specific Items\\label{attrition_covs}",
               coef_map = c(male = "Male",
                            age_num = "Age",
                            religiosity_num = "Religiosity",
                            therm_arabs_num = "Arab Therm"),
               gof_map = c("nobs", "r.squared"),
               stars = c("*"=0.05),
               notes = c("We focus on respondnets who reported covariates in the baseline survey."))  %>%
  add_header_above(c("Outcome: Non Response = 1, Response = 0" = 6))


##Table S4-----
therm_result_age <-
  lm_robust(stdize(therm_arabs_post_num) ~
              pooled*age_num,
            fixed_effects = block_id,
            data = stdy1_dta)

empathy_result_age <- 
  lm_robust(stdize(empathy_index_post) ~
              pooled*age_num,
            fixed_effects = block_id,
            data = stdy1_dta)

peace_result_age <-
  lm_robust(stdize(arabs_peace) ~
              pooled*age_num,
            fixed_effects = block_id,
            data = stdy1_dta) 

exposure_attitudes_result_age <-
  lm_robust(stdize(exposure_index_post) ~
              pooled*age_num,
            fixed_effects = block_id,
            data = stdy1_dta) 

exposure_bhvr_result_age <-
  lm_robust(stdize(bhvr_indx) ~
              pooled*age_num,
            fixed_effects = block_id,
            data = stdy1_dta) 
hte_age <- 
  list("Therm" = therm_result_age, 
       "Empathy" = empathy_result_age, 
       "Peace" = peace_result_age,
       "Support Exposure" = exposure_attitudes_result_age,
       "Exposure Behavior" = exposure_bhvr_result_age)

modelsummary(hte_age,
               title = "Moderating Effect of Age on Primary Outcomes \\label{hte_age}",
               coef_map = c("pooled" = "Treatment",
                            "age_num" = "Age",
                            "pooled:age_num" = "Treat*Age"),
               gof_map = c("nobs", "r.squared"),
               stars = c("*" = 0.05))


##Gender----
therm_result_gnd <-
  lm_robust(stdize(therm_arabs_post_num) ~
              pooled*male,
            fixed_effects = block_id,
            data = stdy1_dta)

empathy_result_gnd <- 
  lm_robust(stdize(empathy_index_post) ~
              pooled*male,
            fixed_effects = block_id,
            data = stdy1_dta)

peace_result_gnd <-
  lm_robust(stdize(arabs_peace) ~
              pooled*male,
            fixed_effects = block_id,
            data = stdy1_dta) 

exposure_attitudes_result_gnd <-
  lm_robust(stdize(exposure_index_post) ~
              pooled*male,
            fixed_effects = block_id,
            data = stdy1_dta) 

exposure_bhvr_result_gnd <-
  lm_robust(stdize(bhvr_indx) ~
              pooled*male,
            fixed_effects = block_id,
            data = stdy1_dta) 



hte_gnd <- 
  list("Therm" = therm_result_gnd, 
       "Empathy" = empathy_result_gnd, 
       "Peace" = peace_result_gnd,
       "Support Exposure" = exposure_attitudes_result_gnd,
       "Exposure Behavior" = exposure_bhvr_result_gnd)

modelsummary(hte_gnd,
               title = "Moderating Effect of Gender on Primary Outcomes \\label{hte_gender}",
               coef_map = c("pooled" = "Treatment",
                            "male" = "Male",
                            "pooled:male" = "Treat*Male"),
               gof_map = c("nobs", "r.squared"),
               stars = c("*" = 0.05))


##Figure S4----
### Personal vs. Control------
therm_result_per_cnt <-
  lm_robust(stdize(therm_arabs_post_num) ~
              personal_vs_control + arab_therm_pre + arab_therm_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Feeling \nThermometer",
         group = "Attitudes",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")


empathy_result_per_cnt <- 
  lm_robust(stdize(empathy_index_post) ~
              personal_vs_control + empathy_index_pre + empathy_index_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Empathy \nIndex",
         group = "Attitudes",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")

peace_result_per_cnt <-
  lm_robust(stdize(arabs_peace) ~
              personal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Peace",
         group = "Attitudes",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")

exposure_attitudes_result_per_cnt <-
  lm_robust(stdize(exposure_index_post) ~
              personal_vs_control + exposure_positive_pre +exposure_positive_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Exposure \nIndex",
         group = "Attitudes",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")

video_result_per_cnt <-
  lm_robust(stdize(video) ~
              personal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Partner \nOrganization \nVideo",
         group = "Behaviors",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")

follow_result_per_cnt <-
  lm_robust(stdize(page_like) ~
              personal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Follow \nPartner \nOrganization",
         group = "Behaviors",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")

exposure_bhvr_result_per_cnt <-
  lm_robust(stdize(bhvr_indx) ~
              personal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Exposure \nBehavior \nCombined",
         group = "Behaviors",
         treatment = "Personal vs. \nControl") %>% 
  filter(.,
         term == "personal_vs_control")



### Non-Personal vs. Control------
therm_result_com_cnt <-
  lm_robust(stdize(therm_arabs_post_num) ~
              communal_vs_control + arab_therm_pre + arab_therm_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Feeling \nThermometer",
         group = "Attitudes",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")

empathy_result_com_cnt <- 
  lm_robust(stdize(empathy_index_post) ~
              communal_vs_control + empathy_index_pre + empathy_index_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Empathy \nIndex",
         group = "Attitudes",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")

peace_result_com_cnt <-
  lm_robust(stdize(arabs_peace) ~
              communal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Peace",
         group = "Attitudes",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")

exposure_attitudes_result_com_cnt <-
  lm_robust(stdize(exposure_index_post) ~
              communal_vs_control + exposure_positive_pre +exposure_positive_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Exposure \nIndex",
         group = "Attitudes",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")

video_result_com_cnt <-
  lm_robust(stdize(video) ~
              communal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Partner \nOrganization \nVideo",
         group = "Behaviors",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")

follow_result_com_cnt <-
  lm_robust(stdize(page_like) ~
              communal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Follow \nPartner \nOrganization",
         group = "Behaviors",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")


exposure_bhvr_result_com_cnt <-
  lm_robust(stdize(bhvr_indx) ~
              communal_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Exposure \nBehavior \nCombined",
         group = "Behaviors",
         treatment = "Non-Personal vs. \nControl") %>% 
  filter(.,
         term == "communal_vs_control")




### Organic vs. Control------
therm_result_org_cnt <-
  lm_robust(stdize(therm_arabs_post_num) ~
              organic_vs_control + arab_therm_pre + arab_therm_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Feeling \nThermometer",
         group = "Attitudes",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")

empathy_result_org_cnt <- 
  lm_robust(stdize(empathy_index_post) ~
              organic_vs_control + empathy_index_pre + empathy_index_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Empathy \nIndex",
         group = "Attitudes",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")

peace_result_org_cnt <-
  lm_robust(stdize(arabs_peace) ~
              organic_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Peace",
         group = "Attitudes",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")

exposure_attitudes_result_org_cnt <-
  lm_robust(stdize(exposure_index_post) ~
              organic_vs_control + exposure_positive_pre +exposure_positive_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Exposure \nIndex",
         group = "Attitudes",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")

video_result_org_cnt <-
  lm_robust(stdize(video) ~
              organic_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Partner \nOrganization \nVideo",
         group = "Behaviors",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")

follow_result_org_cnt <-
  lm_robust(stdize(page_like) ~
              organic_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Follow \nPartner \nOrganization",
         group = "Behaviors",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")

exposure_bhvr_result_org_cnt <-
  lm_robust(stdize(bhvr_indx) ~
              organic_vs_control,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Exposure \nBehavior \nCombined",
         group = "Behaviors",
         treatment = "Organic vs. \nControl") %>% 
  filter(.,
         term == "organic_vs_control")


### Personal vs. Non-Personal------
therm_result_per_com <-
  lm_robust(stdize(therm_arabs_post_num) ~
              personal_vs_communal + arab_therm_pre + arab_therm_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Feeling \nThermometer",
         group = "Attitudes",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")

empathy_result_per_com <- 
  lm_robust(stdize(empathy_index_post) ~
              personal_vs_communal + empathy_index_pre + empathy_index_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Empathy \nIndex",
         group = "Attitudes",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")

peace_result_per_com <-
  lm_robust(stdize(arabs_peace) ~
              personal_vs_communal,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Peace",
         group = "Attitudes",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")

exposure_attitudes_result_per_com <-
  lm_robust(stdize(exposure_index_post) ~
              personal_vs_communal + exposure_positive_pre +exposure_positive_pre_na,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Exposure \nIndex",
         group = "Attitudes",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")

video_result_per_com <-
  lm_robust(stdize(video) ~
              personal_vs_communal,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Partner \nOrganization \nVideo",
         group = "Behaviors",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")

follow_result_per_com <-
  lm_robust(stdize(page_like) ~
              personal_vs_communal,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = "Follow \nPartner \nOrganization",
         group = "Behaviors",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")

exposure_bhvr_result_per_com <-
  lm_robust(stdize(bhvr_indx) ~
              personal_vs_communal,
            fixed_effects = block_id,
            data = stdy1_dta) %>% tidy() %>% 
  mutate(.,
         outcome = " Exposure \nBehavior \nCombined",
         group = "Behaviors",
         treatment = "Personal vs. \nNon-Personal") %>% 
  filter(.,
         term == "personal_vs_communal")


disagg_treat <-
  rbind(therm_result_per_cnt, empathy_result_per_cnt, peace_result_per_cnt,
        exposure_attitudes_result_per_cnt, video_result_per_cnt, follow_result_per_cnt,
        exposure_bhvr_result_per_cnt,therm_result_com_cnt, empathy_result_com_cnt,
        peace_result_com_cnt,exposure_attitudes_result_com_cnt, video_result_com_cnt,
        follow_result_com_cnt ,exposure_bhvr_result_com_cnt, therm_result_org_cnt,
        empathy_result_org_cnt, peace_result_org_cnt, exposure_attitudes_result_org_cnt,
        video_result_org_cnt,follow_result_org_cnt, exposure_bhvr_result_org_cnt,
        therm_result_per_com, empathy_result_per_com, peace_result_per_com,
        exposure_attitudes_result_per_com, video_result_per_com, follow_result_per_com,
        exposure_bhvr_result_per_com)


##Plot results------
ggplot(disagg_treat, aes(x=estimate, y=outcome, color = treatment, shape = treatment)) + 
  geom_point(position = position_dodge(width = 0.6))+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1, position = position_dodge(width = 0.6)) +
  facet_wrap(~group,
             scales = "free") +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  scale_color_manual(values = c("dodgerblue4", "firebrick3", "gray30", "seagreen"))+
  labs(x = "",
       y = "",
       color = "",
       shape = "")+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 16),
        # axis.text.y = element_text(size=14),
        legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank()) 



##Table S8----
# Select the variables of interest
covs_fb <- stdy1_dta %>%
  dplyr::select(age, gender, religiosity)

# Calculate proportions excluding "No-Response"
proportions_df <- covs_fb %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value") %>%
  filter(Value != "No-Response") %>%  # Exclude "No-Response"
  group_by(Variable, Value) %>%
  summarise(
    Count = n(),  # Count of each value
    .groups = "drop_last"  # Retain grouping by Variable for next step
  ) %>%
  mutate(
    Total = sum(Count),  # Total counts for each variable
    Proportion = Count / Total  # Proportion for each value
  ) %>%
  ungroup() %>%  # Remove all grouping
  arrange(Variable, desc(Proportion)) %>% 
  dplyr:: select(Variable,Value,Proportion) %>% 
  mutate(.,
         Variable = case_when(
           Variable == "age" ~ "Age",
           Variable == "religiosity" ~ "Religiosity",
           Variable == "gender" ~ "Gender"
         ),
         Proportion = round(Proportion,3))

# Save the LaTeX table
proportions_df %>%
  kable(format = "latex", 
        booktabs = TRUE, 
        label = "tab:desc_1",
        align = "c", 
        caption = "Descriptive Statistics Facebook Experiment")

